This course will be an introduction to Astropy, a maturing library for astronomy routines and tools in Python.
Astropy started as a combination of various common Python libraries (Pyfits, PyWCS, asciitables, and others) and is working towards providing a consistent API with capabilities for all astronomers. It is developed with extensive automated testing, long-term stable releases, extensive documentation, and a friendly community for contributions.
Note that this design differs from the IDL Astronomy User's Library, which is essentially a mishmash of routines.
These tutorials make some use of the examples at:
In [1]:
# First, make sure this works:
import astropy
# If this doesn't work, raise your hand!
In [2]:
# First we load the fits submodule from astropy:
from astropy.io import fits
# Then we load a fits file (here an image from the Schmidt telescope)
hdu_list = fits.open('http://data.astropy.org/tutorials/FITS-images/HorseHead.fits')
print(hdu_list)
The FITS file contains two header data units, a Primary HDU and an ASCII table HDU (see NASA's Primer) for the different types and limitations.
We can use the hdu_list object like a list to obtain information about each HDU:
In [3]:
print(hdu_list[0].data)
print(type(hdu_list[0].data))
print(hdu_list[0].header['FILTER'])
print(hdu_list[0].shape)
In [4]:
# We can also display the full header to get a better idea of what we are looking at
hdu_list[0].header
Out[4]:
Since this is an image, we could take a look at it with the matplotlib package:
In [5]:
# If using ipython notebook:
%matplotlib inline
# Load matplotlib
import matplotlib.pyplot as plt
# Load colormaps (the default is somewhat ugly)
from matplotlib import cm
# If *not* using ipython notebook:
# plt.ion()
plt.imshow(hdu_list[0].data, cmap=cm.gist_heat)
plt.colorbar()
Out[5]:
We can manipulate the HDU's in any way that we want with the astropy.io.fits submodule:
In [6]:
hdu_list[0].data /= 2
hdu_list[0].header['FAKE'] = 'New Header'
hdu_list[0].header['FILTER'] = 'Changed'
print(hdu_list[0].data)
hdu_list[0].header
Out[6]:
Let's write our new FITS file to our local computer. Running pwd will tell us what directory it is saved to.
In [7]:
%pwd
Out[7]:
In [8]:
hdu_list.writeto('new-horsehead.fits', clobber=True)
While a number of ASCII file readers exist (including numpy.genfromtxt, numpy.loadtxt, and pandas.read_*), Astropy includes readers text file formats commonly used in Astronomy.
These are read as an Astropy Table object, which are convertable to numpy arrays or pandas DataFrames. These can contain unit information and there is work on-going to incoporate uncertainities.
In [9]:
# First we load the ascii submodule:
from astropy.io import ascii
example_csv = ascii.read('http://samplecsvs.s3.amazonaws.com/Sacramentorealestatetransactions.csv')
print(example_csv)
In [10]:
# We can also read Astronomy-specific formats.
# For example, IPAC formatted files
example_ipac = ascii.read('http://exoplanetarchive.ipac.caltech.edu/docs/tblexamples/IPAC_ASCII_one_header.tbl')
print(example_ipac)
Tables support many of the same indexing and slicing operations as numpy arrays, as well as some of the higher-level operations of pandas. See the Astropy tutorial for more examples.
In [11]:
from astropy import units as u
In [12]:
# SI, cgs, and other units are defined in Astropy:
u.m, u.angstrom, u.erg, u.Jy, u.solMass
Out[12]:
In [13]:
# Units all have documentation and attributes
print(u.solMass.names)
print(u.solMass.physical_type)
We can create composite units, such as units of acceleration:
In [14]:
u.m / u.second / u.second
Out[14]:
In [15]:
u.pc / u.attosecond / u.fortnight
Out[15]:
In addition to unit manipulation, Astropy has a concept of Quantities - numbers (or arrays) with units:
In [16]:
print(5*u.erg/u.second)
5*u.erg/u.second
Out[16]:
In [17]:
import numpy as np
my_data = np.array([1,2,3,4,5,6]) * u.Hertz
print(my_data)
In [18]:
# Quantities (and their units) can be combined through algebraic manipulation:
new_data = (6.626e-34 * u.m**2 * u.kg / u.second) * my_data
print(new_data)
Since the computer knows the physical types of each unit, it is able to make conversions between them. Let's use this to simplify my_data. The decompose method will try to use the most basic units, while the .si and .cgs will attempt simple representations with those two bases:
In [19]:
print(new_data.cgs)
print(new_data.si)
print(new_data.decompose())
In [20]:
# We can use the to() method to convert to anything with the same physical_type
print(new_data.unit.physical_type)
print(new_data.to(u.joule))
print(new_data.to(u.eV))
In [21]:
# With the to() method, unit changes are relatively straightforward:
(420*u.parsec).to(u.AU)
Out[21]:
Astropy also includes constants in another submodule, astropy.constants. For example, the average magnitude of the gravitational force of the Earth on the Sun, in SI units, is:
In [22]:
from astropy.constants import M_earth, G, M_sun
(G * M_earth * M_sun / u.AU**2).to(u.N)
Out[22]:
Astropy will even convert units that are not physically compatible, if you are explicit about how to do the conversion. For example, the relationship between wavelength and frequency of light is defined by the choice of the speed of light, allowing the conversion of one to the other:
In [23]:
(450. * u.nm).to(u.GHz, equivalencies=u.spectral())
Out[23]:
A very useful trick is that Astropy will even convert units that require extra information to do so. For example, flux density is usually defined as a density with respect to wavelength or frequency, with the two forms convertable via: $$ \nu f_\nu = \lambda f_\lambda$$ To convert between the different definitions of flux density, we merely need to supply the wavelength or frequency used:
In [24]:
f_lambda = (1e-18 * u.erg / u.cm**2 / u.s / u.angstrom)
print(f_lambda.to(u.Jy, equivalencies=u.equivalencies.spectral_density(1*u.micron)))
print(f_lambda.to(u.Jy, equivalencies=u.equivalencies.spectral_density(299.79*u.THz)))
What about units that coorrespond to locations?
While there does exists u.degree and u.arcsecond, the essential coordinate manipulation is part of the astropy.coordinates submodule. Coordinate conversions, catalog conversions, and more are supported.
In [25]:
# Let's import the main class used, SkyCoord, and create a couple SkyCoord objects:
from astropy.coordinates import SkyCoord
print(SkyCoord(-2*u.deg, 56*u.deg))
print(SkyCoord(1*u.hourangle, 5*u.degree))
print(SkyCoord('2h2m1s 9d9m9s'))
print(SkyCoord('-2.32d', '52.3d', frame='fk4'))
print(SkyCoord.from_name("M101"))
sc = SkyCoord('25d 35d')
In [26]:
# We can retrive the coordinates we used to create these objects:
print(sc.ra)
print(sc.dec)
In [27]:
# We can transform coordinates to different representations (i.e., coordinate systems)
print(sc.transform_to('fk4'))
print(sc.transform_to('galactic'))
In [28]:
# Seperations and position angles are calculatable from SkyCoord objects:
sc2 = SkyCoord('35d 25d', frame='galactic')
print(sc.separation(sc2))
In [29]:
# When we have a world coordinate system (e.g., from a FITS file header), we can convert to and from pixel coordinates:
from astropy.wcs import WCS
w = WCS(hdu_list[0].header)
print(sc.to_pixel(w))
print(SkyCoord.from_pixel(5, 5, w))
In [30]:
import numpy as np
from astropy.modeling import models, fitting
np.random.seed(0)
x = np.linspace(-5., 5., 200)
y = 3 * np.exp(-0.5 * (x - 1.3)**2 / 0.8**2)
y += np.random.normal(0., 0.2, x.shape)
plt.plot(x, y, 'ko')
plt.xlabel('Position')
plt.ylabel('Flux')
Out[30]:
In [31]:
# Fit the data using a Gaussian
model_object = models.Gaussian1D(amplitude=1., mean=0, stddev=1.)
fiter = fitting.LevMarLSQFitter()
g = fiter(model_object, x, y)
plt.plot(x, y, 'ko')
plt.plot(x, g(x), 'r-', lw=2)
plt.xlabel('Position')
plt.ylabel('Flux')
Out[31]:
In [32]:
# Get information about the fitted model:
print(g)
In [33]:
# Load the 9-year WMAP Cosmology and get H_0
from astropy.cosmology import WMAP9 as cosmo
print(cosmo.H(0))
In [34]:
# find the age of the universe at a given redshift:
print(cosmo.age(1))
In [35]:
# Other cosmologies are avaliable
from astropy.cosmology import Planck13 as cosmo
print(cosmo.age(1))
In [36]:
# Build your own cosmology
from astropy.cosmology import FlatLambdaCDM
cosmo = FlatLambdaCDM(H0=70, Om0=0.3, Ob0=0.05)
print(cosmo.age(1))
I don't believe I'll have time to address these, but it is worth pointing out these useful packages:
matplotlib.